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Abstract 


Delay differential equations (DDEs) play an important role in the 
modeling of dynamic processes. Delays arise in contemporary control 
schemes like networked distributed control and can cause deterioration 
of control performance, invalidating both stability and safety properties. 
This induces an interest in DDE especially in the area of modeling 
and verification of embedded control. In this article, we present an 
approach aiming at automatic safety verification of a simple class of 
DDEs against requirements expressed in a linear-time temporal logic. 
As requirements specification language, we exploit metric interval 
temporal logic (MITL) with a continuous-time semantics evaluating 
signals over metric spaces. We employ an over-approximation method 
based on interval Taylor series to enclose the solution of the DDE 
and thereby reduce the continuous-time verification problem for MITL 
formulae to a discrete-time problem over sequences of Taylor coefficients. 
We encode sufficient conditions for satisfaction as SMT formulae over 
polynomial arithmetic and use the iSAT3 SMT solver in its bounded 
model-checking mode for discharging the resulting proof obligations, 
thus proving satisfaction of time-bounded MITL specifications by the 
trajectories induced by a DDE. In contrast to our preliminary work in 
[44], we can verify arbitrary time-bounded MITL formulae, including 
nesting of modalities, rather than just invariance properties. 
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1 Introduction 
“Indecision and delays are the parents of failure.” 


[attributed to George Canning, 1770-1827] 


Ordinary differential equations (ODEs) are traditionally used to model the 
continuous behavior within continuous- or hybrid-state feedback control 
systems. Significant research has consequently been pursued to achieve auto- 
matic verification for such dynamical systems, among it seamless integration 
of safe numeric ODE solving with satisfiability-modulo-theory solving [9, 14]. 
In practice, delay is introduced into the feedback loop if components are 
spatially or logically distributed. Such delays may significantly alter the 
system dynamics and unmodeled delays in a control loop consequently have 
the potential to invalidate any stability and safety certificate obtained on 
the delay-free model. An appropriate generalization of ODE able to model 
delays within the framework of differential equations is provided by delay 
differential equations (DDEs), as suggested by [3]. 

DDEs play an important role in the modeling of natural or artificial 
processes with time delays in biology, physics, economics, engineering, etc. 
As a consequence, attention has gone to developing tools permitting their 
mechanical analysis. However, such tools still are mostly confined to numeric 
simulation, e.g. by Matlab’s dde23 algorithm. Numerical simulation, despite 
being extremely useful in system analysis, fails to present reliable certifi- 
cates of system properties due to numeric approximation. Techniques for 
safely enclosing set-based initial value problems of ODEs, be it safe interval 
enclosures [30, 40, 25], Taylor models [4, 31], or flow-pipe approximations 
based on polyhedra [7], zonotopes [15], ellipsoids [22], or support functions 
[24], consequently need to be lifted to DDEs. 

A safe enclosure method using Taylor series with coefficients in interval 
form was presented in [44]. To avoid dimension explosion incurred by the 
ever-growing degree of the Taylor series along the time axis, the method 
depends on fixing the degree for the Taylor series and moving higher-degree 
terms into the parametric uncertainty permitted by the interval form of 
the Taylor coefficients. By using this data structure to iterate bounded 
degree Taylor over-approximations of the time-wise segments of the solution 
to a DDE, the approach identifies the operator that yields the parameters 
of the Taylor over-approximation for the next temporal segment from the 
current one. Employing constraint solving to analyze the properties of this 
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operator, an automatic procedure is obtained to provide stability and safety 
verification for a simple class of DDEs of the form 


< a(t) = f(#t-8)) () 

with linear or polynomial vector field f : RN —» R‘, where the derivative 
at t is a function of the trajectory at t — 6, i.e., the signal value determines 
the future evolution with delay of 6. The limitation of the method proposed 
in [44] is that its coverage of safety properties is confined to the verification of 
state invariants only. Improving on the previous work in [44], the contribution 
of the current article lies in verifying a class of safety requirements specified 
using linear-time temporal logic. 

The method proposed in this article again addresses DDEs in the form 
of Eq. (1) and builds upon the safe enclosure method for DDEs presented in 
[44], yet addresses verification of behavioral properties expressed in metric 
interval temporal logic (MITL) [1, 11, 33]. MITL is a linear temporal logic 
that is meaningful when the states evolve in metric spaces, an assumption 
met by continuous-state systems as in Eq. (1). It is considered as a real- 
time extension of linear temporal logic (LTL), where the modalities of LTL 
are constrained with timing bounds. In particular, given a continuous 
dynamical system (1) with its initial condition(s) and a temporal logic 
specification expressed in time-constrained MITL, we employ the interval- 
based Taylor over-approximation method to enclose the solution of the given 
DDE. This facilitates effective reduction of the signal-based, continuous-time 
and continuous-state MITL verification problem to a related discrete-time 
MITL verification problem expressible in terms of timed state sequences. 
By using any bounded model checking (BMC) tool built on top of an 
arithmetic SMT solver being able to address polynomial arithmetic, we obtain 
a procedure able to provide safety certificates for DDE relative to temporal 
logic specifications. In our case, we use the iSAT3? implementation of the 
iSAT algorithm [13] that provides techniques for bounded and unbounded 
verification problems like k-induction [39] and Craig interpolation [29]. 

For dealing with temporal properties expressed in MITL, the key step 
is to safely determine truth values of atomic propositions, i.e., to generate 
sufficient conditions for their validity over a time frame based on the Taylor 
over-approximation of the DDE (1). Based on this, the solver is able to verify 
more complex formulae of temporal logic also involving Boolean connectives 
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and temporal modalities, like the (bounded) until operator. Our approach 
is characterized by the soundness guarantees obtained due to the over- 
approximation of the DDE and the sufficient conditions used for substituting 
the exact MITL semantics. The accuracy of approximation can be selected; 
an automatic refinement method dynamically adapting the accuracy in case 
of a negative verdict, however, remains to be developed. We demonstrate 
how our approach works in practice by presenting verification of temporal 
properties on example systems. 


Structure of the article. After discussing related work in Section 2, 
we formulate the temporal verification problem on DDE in the form of 
Eq. (1) by defining syntax and continuous-time, continuous-state signal- 
based semantics of MITL formulae, our requirements specification language 
(Sect. 3). Section 4 develops interval-based Taylor over-approximation as a 
safe time-wise discretization to the solution of the DDEs, providing a time- 
invariant operator generating a timed state sequence on Taylor coefficients. 
In Section 5, we adapt the interpretation of MITL to the timed state sequence 
such that it safely recovers the original semantics on the actual solution of 
the DDE in terms of conditions on the Taylor coefficients of the time-discrete 
model. Finally, we conclude our paper in Section (6) presenting some ideas 
for further refinement and future directions of our work. 


2 Related Work 


Driven by the demand for safety cases (in a broad sense) for safety-critical 
control systems, we have over the past decades seen a rapidly growing 
interest in automatic verification procedures for system models involving 
continuous quantities and dynamics described by, a.o., differential equations. 
Traditionally ordinary differential equations (ODEs) are used for describing 
system dynamics and safety properties have been specified in terms of a set 
of unsafe states formulating the verification problem as reachability analysis. 
Reachability analysis, which involves computing appropriate approximations 
of the reachable state sets, plays a fundamental role in addressing such safety 
verification challenges. Consequently, significant research has been invested 
in reachability analysis of such dynamical systems [30, 4, 31, 7, 15, 22, 24]. 
In contrast to this extensive literature on computing the reachable state 
space for ODEs, the reachability analysis to dynamic systems modelled by 
delay differential equations (DDEs) is in its infancy and thus provides an 
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open area of research. 

Zou, Franzle et al. proposed in [44] a safe enclosure method using 
interval-based Taylor over-approximation to enclose a set of functions by a 
parametric Taylor series with parameters in interval form for a simple class of 
DDEs in the form of Eq. (1). This method dealt only with simple invariants 
as safety properties. In [35], Prajna et al. extended the barrier certificate 
methodology for ODEs to the polynomial time-delay differential equations 
setting, in which the safety verification problem is formulated as a problem 
of solving sum-of-square programs. The work in [18] presents a technique for 
simulation-based time-bounded invariant verification of nonlinear networked 
dynamical systems with delayed interconnections by computing bounds on 
the sensitivity of trajectories (or solutions) to changes in initial states and 
inputs of the system. A similar simulation method integrating error analysis 
of the numeric solving and the sensitivity-related state bloating algorithms 
was proposed in [6] to obtain safe enclosures of time-bounded reach sets for 
systems modelled by DDEs. 

Confining safety properties to a set of unsafe states, as in the afore- 
mentioned work, considerably restricts the ability of designers to adequately 
express the desired safe behavior of the system that may involve a number 
of critical properties such as timing requirements and bounded response. 
Metric temporal logic (MTL), introduced by Koymans [20], is popular for- 
malism for expressing such properties as a real-time extension of linear 
temporal logic (LTL) [28] to specify real-time properties. Then, Alur et 
al. in [1] introduced metric interval temporal logic (MITL) to address the 
undecidability problem of MTL by relaxing the punctuality of the temporal 
operators. The bounded-time verification or falsification of such properties 
has been studied for continuous/hybrid systems in [10, 11, 37, 34, 27], yet 
DDEs are not handled. In this paper, however, we present an approach 
aiming at automatic safety verification of a class of DDEs in the form of 
Eq. (1) against requirements expressed in MITL formulae. In contrast to 
our preliminary work in [44], yet we can verify arbitrary time-bounded 
MITL formulae including nesting of modalities rather than just invariance 
properties. 


3 Problem Formulation 


In this section, we formulate the verification problem of a simple class of 
DDEs in the form of Eq. (1) against a class of safety requirements specified 
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using an appropriate linear-time temporal logic. As we deal with continuous 
state and time, we adopt metric interval temporal logic (MITL) [1] for the 
purpose of requirements specification language. In this section, we review 
its syntax and its continuous-time, signal-based semantics. 

Let R be the set of the real numbers. Our time domain is the set 
of nonnegative real numbers Rs. Also, the trajectory of the DDE of Eq. 
(1) on an initial condition x({0,6]) = c € R is a function x(t) such that 
az: Ryo > RN satisfies the initial condition and Vt > 6: {2z(t) = f(Z(t—64)), 
where the positive integer N denotes the dimension of the state space. In 
order to specify the temporal properties of interest, we exploit MITL with 
continuous semantics, as meaningful when the states evolve in metric spaces 
like in Eq. (1). We say that P(C) denotes the powerset of a set C and assume 
that AP is a set of atomic propositions. Then, the predicate mapping 
M : AP + P(RY) is a set valued function that assigns to each atomic 
proposition p € AP a set of states M(p) C R™. In this paper, we take 
the set of atomic propositions AP to be bound constraints e ~ c on state 
expressions, where e is an expression formed over the state variables, like 
1x2 — 2sinx3, and being compared via a relation ~ € {<,<,>,>} toa 
constant c € Q. Such atomic propositions come equipped with their natural 
semantics. 


3.1 Metric Interval Temporal Logic 


Metric interval temporal logic (MITL) [1] is a linear-time temporal logic 
designed for capturing properties of signals evolving over quantitative and 
thus metric rather than qualitative time, an assumption met by continuous- 
state systems as in Eq. (1). It is a real-time extension of linear temporal logic 
(LTL), where the modalities of LTL are constrained with quantitative timing 
bounds. Metric temporal logic (MTL) was first introduced by Koymans 
[20] to specify real-time properties. In order to address the undecidability 
problem of MTL, Alur et al. in [1] relaxed the punctuality of the temporal 
operators s.t. they cannot constrain to singleton intervals. We employ MITL 
to formally characterize the desired behavior of DDEs. Along the following 
lines, we review and suitably adapt the syntax and the continuous-time and 
continuous-space semantics of MITL as presented in [1, 33]. 


Definition 1 (Syntax of MITL). An MITL formula ¢ is built from a set 
of atomic propositions AP using Boolean connectives and timed-constrained 
versions of the until operator. It is inductively defined according to the 
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grammar 


gu=T |p| 7g. | ¢1A 2 | ¢1 Uz 2 


where p € AP, T is the Boolean constant true and Z C Q>so is a nonsingular 
interval imposing timing bounds on the temporal operators, where Q>o is the 
set of non-negative rational numbers. 


We can derive the constant false by L = =T. Also, we can define 
additional, time-constrained version of, temporal operators such as release 
Rr, eventually Oz, and always Oz as follows: 


yi Rr pe = 7((741) Uz (-¢%2)), 
Or p=TUrz y, and 


TP= Rr p= 707 7. 


Notice that the release operator is a temporal modality that is dual 
to the until operator. A formula y; Rz v2 holds if yo always holds, a 
requirement that is released as soon as y, becomes valid w.r.t. the time 
bounds 7. 

Also, note that MITL has no nezt operator as the time domain is dense. 
For Z = [0, co], we can remove the subscript Z from the temporal operators, 
obtaining the traditional modalities of LTL. Finally, we would like to point 
out that the decidability problem of MITL in the continuous semantics for 
both model checking and satisfiability problems is out of the scope of this 
paper. For details about the decidability problem, we refer the reader to 
[1, 32]. Furthermore, it is an open issue whether the model property of DDE 
w.r.t. MITL formulae is decidable. 


3.1.1 Negation Normal Form 


We consider MITL formulae in negation normal form (NNF'}), which can be 
achieved by pushing all negations inside into the atoms [36]. If we admit the 
release modality and disjunction in our syntax, then every formula y has a 
semantically equivalent negation normal form nnf (py). Such an NNF can be 
obtained by applying De Morgan’s laws as well as the dualities between until 
and release in order to push negations inwards, and thereafter eliminating 
double negations. This is done by exploiting the following equivalences as 
rewrite rules from left to right: 
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7771 = $1; 
a(p1 A g2) = 791 V mH2, 
a(p1 V 2) = 791 A mH, 
a(~1 Ur p2) = 791 Rr 72, 
(yi Rr pe) = 71 Ur ~~. 


These rewrite rules can also be lifted to the derived operators as follows: 


Or p = Ur 79, 
= Tp=Or 7H. 


3.1.2. Continuous-Time, Continuous-State Semantics of MITL 


The continuous semantics of MITL formulae is used to express specifications 
on the desired temporal evolution to the solutions of DDEs in the form of 
Eq. (1). This semantics is based on real-valued signals x : Ryo > R over 
time. We say that expression e over the state variables x satisfies atomic 
formula e ~ c at time t > 0, denoted e,t E e ~ ¢, iff e(t) ~ c holds. Based 
on this, semantics of arbitrary MITL formulae is defined inductively, with 
the semantics of Boolean connectives = and A as well as the constant T 
being standard. The semantics of the time-constrained until operator is 
defined as follows: e,t — yi Uz 2 iff for some t’ € T, e,t +t’ — ve holds 
and furthermore e,t | y; for allt € (t,t +t’). 

By convention, we say that the DDE of Eq. (1) with an initial value 
x({0,6]) = c satisfies an MITL formula if the expression e(t) over its 
solution trajectory satisfies y in the sense of e,0 E y. In what follows, we 
employ the interval-based Taylor over-approximation method from [44] to 
enclose the solution of such a DDE. As this method factually generates a 
discrete sequence of Taylor coefficients rather than a continuous trajectory, 
we are thus able to reduce a correctness problem over continuous time into 
a corresponding problem of a time-invariant operator over discrete time. 
Therefore, it is however necessary to recover the continuous semantics on 
the actual solution of the DDE from the timed state sequence semantics on 
the Taylor coefficients. 


Model Checking Delay Differential Equations 
Against Metric Interval Temporal Logic 85 


4 Computing Enclosures for DDEs by Taylor Mo- 
dels 


In this section, we review the bounded degree interval-based Taylor over- 
approximation method for a simple class of DDEs first presented in [44]°. In 
order to compute an enclosure for the trajectory x(t) defined by an initial 
value problem of the DDE (1), a template interval Taylor form of fixed 
degree & is defined as 


falt) = Gn + Gn, t+---+Gn,t*, (2) 


where f,, encloses the trajectory for time interval [nd, (n + 1)d], the constant 
6 is the feedback delay from Eq. (1), and ap,,...,@n, are interval-vector 
parameters. The trajectory induced by DDE (1) can be represented by 
a piece-wise function, with the duration of each piece being the feedback 
delay 6. To compute the enclosure for the whole solution of the DDE, we 
need to calculate the relation between the interval Taylor coefficients in 
successive time steps as pre-post-constraints on these interval parameters. For 


notational convenience, we denote the interval parameters [an,.,...,@n,] by 
a matrix A(n) in RNx(k+)) The relation between A(n) and A(n +1) can be 
computed, exploiting different orders of Lie derivatives eae age 25 ve 
as follows: 
af), (t) af) 
1 2 n k n 
frrilt) =9(falt))s fen) =... OQ =—AL, (8) 


i.e., the first order is obtained directly from the given DDE (1) and the 
(i + 1)-st order is computed from the i-th order by symbolic differentiation. 
Then, the Taylor expansion of fn+1(t) with fixed degree k is derived as 
follows: 


AGS 
(k — 1)! "kl 


ft (0) fra ©) 4 


fn+i(t) — fn (9) a 7 bap eseae ae (4) 


where €, is a vector ranging over [0, 6]. 


3The corresponding prototype implementation of the interval Taylor over-approximation 
method for DDEs as well as some examples are available for download from 
https: //github.com/liangdzou/isat-dde. 
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From Eq. (4), by comparing the coefficients of monomials with the same 
degree at the two sides and by replacing €, by the interval vector [0,6], we 
can obtain a time-invariant operator which represents the relation between 
A(n) and A(n +1). The details of this construction can be found in [44] 
or retrieved from the example underneath. Hence, we safely enclose the 
trajectory induced by the DDE (1) by a discrete-time model providing a 


timed state sequence on a state space S C RN*(#+)), 


4.1 Time-Wise Discretization of DDEs into Timed State Se- 
quences 


We demonstrate on a running example taken from [44] how to provide the 
discrete-time model that encloses the solution of a DDE like Eq. (1). The 
running example is the DDE 


#(t) = —2(t - 1) (5) 


with the initial condition x([0,1]) = 1. Fig. 1 shows the solution of ODE 
x = —x without delay (the dashed line) and with 1 second delay (the solid 
line). Obviously, the difference between the ODE and the DDE is substantial 
and necessitates analysing the behaviour of the DDE. 


1 


0.54 


-0.5- 


Figure 1: Solutions to the ODE « = —2 (dashed graph) and the related DDE «(t) = 
—ax(t — 1) (solid line), both on similar initial conditions (0) = 1 and 2([0,1]) = 1, 
respectively. 
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The method provided in [44] aims at over-approximating the solu- 
tion of DDE (5) by iterating bounded degree interval-based Taylor over- 
approximations of the time-wise segments of the solution to the DDE. That 
way, we identify the operator that yields the parameters of the Taylor over- 
approximation for the next temporal segment from the current one. For 
instance, suppose we are trying to over-approximate the solution of DDE (5) 
by polynomials of degree 2. Then we can predefine a template Taylor form 
fn(t) = no + Qn,t + Qnyt? on interval [n,n +1], where an, @n,, and Gn, are 
interval parameters able to incorporate the approximation error eventually 
necessarily incurred by bounding the degree of the polynomial to (in this 
example) 2. Here, f,,(t) corresponds to the solution x of DDE (5) at time 
n+t, ie., f(t) over-approximates x(n +t) in the sense of x(n + t) € f(t). 
In order to compute the Taylor model, the first and second derivative 


» (t) and f(t) of solution segment n+1 based on the preceding segment 


(where both segments are of duration 1 each) have to be calculated. The 


a 


first derivative FO (0) is computed directly from Eq. (5) as 
fret (t) = —Sn() = Ong — Ont — Ange? 


The second derivative f(t) is computed based on fi) by 


2 
F208) = = —An, — 2anyt. 
By using a Lagrange remainder with fresh variable €, € [0,1], we obtain 


fues(0), . Li 


2 
1! 2! : 


Fati(t) = fn(1) 4 


any + 2an2En t2 
5 : 


= (ng + An, + Anz) — Anot 


Then, the operator expressing the relation between Taylor coefficients 
in the current and the next step can be derived by replacing both f;,(t) and 
fn+i(t) with their parametric forms an. + nit + @n.t? and @n419 + Qn41,t + 
Gn+1,t” in the above equation and pursuing coefficient matching. As a result, 
one obtains the operator 


An+lo 1 1 1 Ang 
Qn41i,| = |-1 0 0 | Jan, (6) 


An+lo 
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mapping the coefficients of the Taylor form at step fn to the coefficients of 
the Taylor form of f,41. The coefficients change at every 6 time units (every 
second in the given example) according to the above operator, which therefore 
defines a discrete-time dynamical system corresponding to the DDE. The 
discrete-time operator can be rendered time-invariant, yet interval-valued 
by substituting the uncertain time varying parameter €,, with its interval 
[0,6]. Hence, we can safely enclose the solution of DDE (5) by a sequence 
of parametric Taylor series with parameters in interval form. In the case of 
system (5), as well as for any other linear DDE, the operator generating this 
sequence is a set-valued linear operator definable by an effectively computable 
interval matrix. 


4.2 Bounded Model Checking Mode in iSAT3 


In order to encode the Taylor model corresponding to a DDE, we use 
bounded model checking (BMC) mode in iSAT3 [38]. The iSAT3 solver is 
a satisfiability checker for Boolean combinations of arithmetic constraints 
over real- and integer-valued variables as well as a bounded model-checker 
for transition systems over the same fragment of arithmetic. It is a stable 
version implementation of the iSAT algorithm [13]. The solver can efficiently 
solve bounded verification problems that involve polynomial (and, if needed, 
transcendental) arithmetic. Hence, it is a good option to solve our proposed 
problem due to the Taylor forms involved. Also, it allows us to verify/falsify 
a variety of MITL formulae built on atomic predicates defined over simple 
bounds, linear, and nonlinear constraints [21]. Bounded model checking 
(BMC) of a transition system aims at finding a run of bounded length kgeptn 
which 


e starts in an initial state of the system, 
e complies with the system’s transition relation, and 
e ends in a state in which a certain (un)desired property holds. 


The bounded model checking engine then constructs a formula which is 
satisfiable if and only if a trace with above properties exists.4 In case of 
satisfiability, any satisfying valuation of this formula corresponds to such a 


‘Tt should be noted that this semantic property does not imply that the solver engine 
subsequently checking that formula for satisfiability can exactly determine its satisfiability. 
In the case of iSAT, a sound, yet incomplete unsatisfiability check is implemented, as 
necessitated by the undecidable fragment of arithmetic addressed. 
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trace. For encoding the discrete transition system on Taylor model in BMC 
mode, iSAT3 has an input file format consisting of four sections: 


e DECL: This part contains declaration of all variables (i.e., variables of 
the dynamic system, Taylor coefficients of the Taylor over-approximation 
solution, the duration of each segment t € [0,6], the uncertain time- 
varying parameter € € [0,6]). 


e INIT: This part is a formula describing the initial state(s) of the system 
to be investigated. 


e TRANS: This formula describes the transition relation in symbolic form; 
in our case the evolution of the time-discrete Taylor model. We encode 
a template interval Taylor form of fixed degree k, i.e., fn(t), and the 
relation between interval Taylor coefficients in the current and the next 
step. Variables may occur in primed (e.g., a’) or unprimed (e.g., a) 
form. A primed variable represents the value of that variable in the 
successor step, i.e., after the transition has taken place. 


e TARGET: This formula characterises the state(s) whose reachability is 
to be checked; in our case it represents satisfaction of the given MITL 
formula. 


The solver unwinds the transition relation kgeptn times, conjoins the resulting 
formula with the formulae describing the initial state(s) and the target 
state(s), and then solves the obtained formula. For our transition relation 
in terms of Taylor coefficients, the solver recursively for each time frame 
[0, kaeptnd] constructs the following formula: 


kdepth—1 
init (a°) A \ trans (a, a’*!) A target (Aker), 
i=0 


where a is the interval-vector of the Taylor coefficients of the fixed-degree 
Taylor polynomial f,,(t). 


Back to our running example in Sect. 4.1, the iSAT3 encoding for 
this example is as shown in Listing (1). In the DECL part, we declare all 
variables; the variables of the dynamic system, i.e., x, the Taylor coefficients 
of degree 2, i.e., ag,a 1, and ag, the duration of each segment t € [0,1], and 
the uncertain time varying parameter € € [0,1]. Notice that the range of 
each variable has to be bounded in iSAT3. We initialize the system variable 


SAIDWIKWNH 
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x and the Taylor coefficients in INIT part according to the given initial 
condition(s) in our example. Then, in the TRANS part, we state the interval 
Taylor form of degree 2, i.e., f;,(t) corresponds to the solution x of DDE 
(5), as shown in line 22, and the relation between Taylor coefficients in 
the current (unprimed variables) and the next segment (primed variables) 
according to the generated operator (6), where the segments are of duration 
1 each. 


DECL 

-- the range of each variable has to be bounded 
float [-1000, 1000] a0, al, a2, x; 
float [0,1] t, xi; 


INIT 
-- initial value of solution 


x = 1; 


-- initialize Taylor coefficients 


aQO = x; 

ail = 0; 

a2 = 0; 
TRANS 


-- relation betw. Taylor coefficients in current and next step 
a0’? = aO + al + a2; 
ai’ = -a0; 
a2’? = -0.5*al - xix*a2; 


-- x(t) is given by a Taylor form of degree 2 


x? = a0’ + al’*t + a2’*(t72); 
TARGET 
-- state to be reached, e.g. 

x = -0.25; 


Listing 1: The encoding in iSAT3 of the running example in Sect. 4.1. 


4.3. Proving Continuous-Time Properties on the Time Dis- 
cretization 


Operator (6) straightaway defines a safe temporal discretization of the DDE 
system in Eq. (1), i-e., an operator generating a classical timed state sequence 
in the sense of [1, 10]. We can, however, not simply apply the discrete-time 
interpretation of MITL to this timed state sequence, as it ranges over a 


Model Checking Delay Differential Equations 
Against Metric Interval Temporal Logic 91 


different state space —namely the Taylor coefficients— than the requirements 
specification in terms of the original state variables. Therefore, we have to 
translate forth and back between the different state spaces and time models. 
In detail, the iterated execution of operator (6), starting from an initial 
vector do,,---, 40, Of Taylor coefficients encoding the initial solution segment 
x({0,6]), generates a timed state sequence over (interval) Taylor coefficients, 
with time stamps t; = i6, rather than a continuous signal over the state 
variables 71,...,¢y. Reflecting this encoding, we need a translation step 
generating conditions over the timed sequence of Taylor coefficients from 
which we are able to recover the original continuous-time, continuous-state 
signal-based semantics on the actual solution x of the DDE, as defined in 
Sect. 3.1. 

As has already been observed in [44], such a mapping is straightforward 
when invariance properties are to be dealt with, for which a sufficient —yet, 
in the light of over-approximation of the solution, obviously not necessary— 
condition can be obtained as follows. For an invariance requirement Ox € 
Safe, where Safe is a set of safe states, the requirement in the n-th segment 
is translated to the stronger condition 


Vt € [0,6] VE € [0,6] Vao € Ano,...,a% € Ann: fn(t) € Safe, (7) 


where f, is the underlying Taylor form and Ayo to An, are the intervals 
Taylor coefficients stemming from the n-th iteration of the operator (6). As 
this Taylor form provides an over-approximation of the solution x over time 
frame [nd, (n+ 1)d], the condition (7) implies Vt € [nd, (n+1)d] : x(t) € Safe. 
Consequently, the continuous-time safety property Ox € Safe for system (5) 
is translated into a sufficient condition according to Eq. (7) for t,€ € [0, 1] 
over the sequence of Taylor coefficients of Taylor polynomial of degree 2. 
As its violation is an existential statement both instantiations of Taylor 
coefficients within given intervals and existentially quantified time points t 
and €, a solver for satisfiability modulo theory over the existential theory of 
polynomial arithmetic can be used to solve the safety verification problem. 
It requires polynomial constraint solving due to the Taylor forms, i.e., 
polynomial expressions involved in the statement f,,(t) € Safe. 

Different proof schemes can be implemented using such a solver: using k- 
induction [39] or interpolation-based unbounded proof schemes [29], absence 
of any time point in the sequence of valuations generated by operator (6) 
satisfying dn € N, at € [0,1], 5€ € [0,1],54ao € Ano,...,a% € Anz: fn{t) ¢ 
Safe can be shown, thereby rigorously showing safety of the DDE system 
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under investigation. Bounded model checking of the same system could, on 
the other hand, generate counterexamples to safety, which may however be 
spurious due to the over-approximation involved in the Taylor enclosure. 


5 Solving Continuous-Time MITL Formulae by Re- 
duction to Time-Discrete Taylor Approximati- 
ons 


We extend the above idea of generating sufficient conditions for MITL 
specifications on DDEs in terms of the sequences of enclosing (interval) 
Taylor coefficients. The aim is to cover a large fragment of MITL, expanding 
well beyond the invariance properties addressed in [44]. As explained in the 
previous section, we have obtained a generator for a timed state sequence 
—the operator (6)— representing the solution of the DDE, yet ranging 
over a different state space, namely the Taylor coefficients. Hence, the 
continuous interpretation of the MITL formulae over DDE solutions has 
to be translated into a semantically appropriate discrete interpretation on 
a timed state sequence with time stamps t; = id. This translation needs 
to restore, in the sense of providing sufficient conditions for the solution 
being a counterexample (i.e., a witness of violation of the property), the 
continuous semantics of the MITL formulae over the discrete model of the 
timed state sequence. We do so by first transforming the MITL formula into 
negation normal form, then generating a sufficient condition by adding the 
appropriate conditions to the Taylor model that meet the semantics of the 
property for searching for (possibly spurious) counterexamples with the help 
of an SMT solver. 


5.1 Atomic Proposition 


According to the MITL syntax of Sect. 3, atomic propositions are of the form 
e ~ c, where € is an expression over the state variables, c is a constant and 
~ an inequational relational operator, i.e., one of <,<, >, >. Using bounded 
model-checking based on SMT solving, we attempt to find a counterexample 
of the MITL formula, or, in other words, look for a witness for the negation of 
the MITL formula. As we transform that negated formula into NNF, atomic 
propositions occur in positive context only. Then, sufficient conditions for 
truth of such propositions throughout a time frame [id, (¢ + 1)d6] can, as 
already observed in Eq. (7)), obviously be expressed as follows: 
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Vte (0, 6] VEE [0, 6] Vag — Aio, ++, Ok E Aik ; \ Y= fi(t) Newc. (8) 
41 


As mentioned in Sect. 4.3, when using SMT solving for finding violations 
of condition (8), we use the negation of the universally quantified condition 
Eq. (8). As this is an existential formula, it is amenable to standard SMT 
solving. 


5.2 Boolean Connectives 


For solving complex-structured formulae, we use a Tseitin-like definitional 
translation [43], where we introduce a fresh Boolean helper variable (w); 
for each subformula ~ and each index i of a time frame [id, (i + 1)6]. The 
intuition is that (7); being true implies that w holds for each time point 
t € [id,(¢ + 1)d]. Note that this is a one-sided implication, as we cannot 
decide properties exactly. 

Note that we have in the previous section already obtained appropriate 
definitions for the case that ~ is an atomic formula, such that we can define 


t € [0, 6] FE € [0,4] 7 (9) 


os ( Nya & = Filt) Ae x C 


as sufficient condition for validity of an atomic formula e ~ c, where e is an 
expression over the state variables and ~ is the converse of the relation ~. 
Given a compound formula of the form a, = y1 A Yo or wo = 1 V Ya, 
the encoding for the compound formula is obtained by conjoining to the 
“axiomatisations” of y; and v2 the following definitional translations: 


(p1 A po)i = (y1)i A (2): 
(1 V po)i = (y1)i V (Ya): 


Note that a single-sided implication “=” from right to left would actually 
suffice, as we target sufficient conditions only. 


5.38 Unary Temporal Operators 


Assume we have an MITL formula wy, = Oz y or yo = Oz y featuring a 
time-constrained eventually or always temporal operator as its outermost 
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operator. Let the lower and upper bound of Z for simplicity be integer 
multiplies 16 and ué of 6. For each time frame, the value of a given MITL 
formula is encoded with the help of new Boolean variables for the truth 
values of its subformulae in particular time instants. The encoding of 
and w2 can be recursively understood as follows: 


(O116,u5])8 > (Cfo,(u—Ya] Pi fi<lea 
(Ofousy Pre > (P)E V (Pfo,(u—1)5] P41 ifl<u 
(Opo,s1P)i > (Y)i 

(Oi5,u4))i > (Ffo(u—ya]P) itt» fl <l<a 
(Fio,us]¥)i = (i A (Fo, (u—1)5]Y) i+ 1» ifl<u 

(Tho, Pi = (Y)i 


Single-sided implications “=” from right to left would again suffice for a 
sound definitional translation. 

In the case of the eventually modality, the condition for detecting 
satisfaction of y is somewhat stronger than necessary, actually requiring it 
to hold throughout a full time frame rather than just once inside. 


5.4 Binary Temporal Operators 


Assume we have a subformula of shape #1 = y1 Uz v2 or Wo = Y1 Rr ve 
featuring a time-constrained until or release operator as its outermost con- 
nective. For simplicity, we assume that the lower bound of Z is 0. Such 
a form can always be achieved by prepending the modality with a unary 
temporal operator. Then the encoding of a sufficient condition for validity 
of w1 or we, resp., over time frame [id, (2 + 1)d] is as follows: 


(G1 Ufo,us] $2)i > (P2)i V (Hide A (61 Ufo,(u—1)5] P2)it1), if l <u 
(p1 Uo] P2)i > (P2)i 

= (p2) 
) 


) 


i 
)i 
(G1 Rio,us] 92)i > (2)i A ((P1)4 V (Y1 Rpow—1)3) Y2)i41), ifl <u 
(Y1 Rios] P2)i > (2): 


As in the case of the eventually modality, the condition for detecting yo 
in the case of until and of v1, resp., in the case of release again is somewhat 
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stronger than necessary, requiring it to hold throughout the respective time 
frame instead of just once inside. 


5.5 Correctness 


Let ~% be an MITL formula and [wo be the definitional translation of w 
obtained by recursively unfolding and conjoining the above definitions of 
(~)o and all (y); occurring therein. Let “(t+ 6) = f(£(t)) be a DDE with 
initial value £([0,5]) = 7, and let A be the interval matrix obtained from it 
due to Eq. (6). Let k be the highest index j of any Tseitin variable (y) ; 
occurring in [7]o. 


Lemma 1 If ay =i A a Git. = AG; A [wo A7(w)9 is unsatisfiable then # 
satisfies w, where Go =i denotes the appropriate initialisation of the Taylor 
coefficients and & is the exact solution of the DDE. 


Proof: The sequence Gop = iA Ns Gj41 = Ad; of interval Taylor forms 
generates an over-approximation of x. The construction of [w]o is such 
that @ =i A ae Gi41 = AG; A [Wlo E (w)o if all trajectories y enclosed 
by the sequence of interval Taylor forms, and thus also z itself, satisfies w. 
Satisfiability of d =i A Nes Gi41 = AG; A [v]o A 7(w)o consequently is a 
necessary condition for violation of w by a. 


Note that dy =i A My G41 = AG; A [W]o A 7(w)o is a purely existential 
statement and thus amenable to standard SAT-modulo-theory solving by 
removing the explicit existential quantifiers in each instant of Eq. (9) by 
introducing fresh variables. 


5.6 Verification Examples 


In this section, we use the iSAT3 SMT solver to discharge the above proof 
obligations. In order to be able to present the encodings in a compact 
form suitable for manual inspection and for publication in print, we slightly 
deviate from a strict implementation of the above scheme, and instead 
employ the bounded model checking (BMC) mode of iSAT and symbolic 
counter variables as abbreviation mechanisms whenever appropriate in the 
search for witnesses as counterexamples of the MITL formulae. The results 
are, however, the same and the logics behind the encodings is equivalent to 
the point it can be in a BMC encoding. In particular, the method for using 
existential arithmetic constraints as sufficient conditions to determine the 
truth values of propositional (sub)formulae based on the over-approximation 
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model of the DDE and thus recover the continuous semantics of the MITL 
formula on the actual solution x of the DDE from the timed state sequence 
semantics is exactly as in Eq. (9). 

We demonstrate the approach based on illustrative examples of DDEs 
in the form of Eq. (1). In our examples, we first consider the DDE (5) 
presented in Sect. 4.1 with different conjectured MITL formulae to be 
verified. Thereafter, we apply our method to an adaptation of Gustafson’s 
model of nutrient flow in an aquarium (three dimensional example) [17, p. 
589f]. 


Example 1 We consider the linear DDE «(t) = —a(t — 1) with initial 
condition x([0,1}) =1 and the conjectured safety property Ojo19) (x < 1.2). 


The bounded degree interval-based sequence of Taylor forms can be generated 
by the operator Eq. (6)). Adopting degree 2 Taylor forms, we can encode 
this generator in the iSAT3 input language as shown in lines 24—26 of Listing 
2. The encoding is a discrete-time dynamic system over the variables x 
representing (snapshots of) the DDE solution, Taylor coefficients of the 
Taylor over-approximation solution, i.e., @9,a1,, and ag, a time point in each 
segment t € [0,1], and the uncertain time varying parameter € € [0,1]. Also, 
we declare a counter to observe the timing bound on the temporal operator. 


In order to solve the given MITL formula Oj9,19) (v < 1.2) in iSAT3 in 
the sense of trying to construct a counterexample, we 


1. in accordance with Eq. (9) search for a time frame within which z, 
being defined as the image of the Taylor polynomial for some t € [0, 1], 
€ € [0,1] in line 29 of the listing, exceeds 1.2, as encoded by condition 
x > 1.2 in the target (line 37), and 


2. enforce the count of the time frame to be at most 9 (target, line 38), 
as time frame n ranges from time n ton +1. 


For verifying the property at hand, it obviously suffices to check this formula 
up to unwinding depth 9.° Such bounds on unwinding depths can in iSAT3 
be set with the --start-depth and --max-depth command line options. 


In our example, the solver outputs that the system is safe for unwinding 
depth 10, i.e., no state satisfying the target property could be reached within 


°iSAT counts unwindings starting from 0 such that an unwinding of depth 9 yields a 
trace comprising 10 time instants. 
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the relevant depth. This constitutes a rigorous proof that the system actually 
satisfies the MITL formula. 


DECL 

-- the range of each variable has to be bounded 
float [-1000, 1000] a0, al, a2, x; 
float [0,1] t, xi; 


-- define counter for the bounded verification problem 
int [0,9] counter; 


INIT 
-- initial value of x over [0,1] 


-- initialize Taylor coefficients 


aO = 1; 
ail = 0; 
a2 = 0; 


-- initialize the counter observing the time interval 
-- covered by the bounded always 
counter = 0; 


TRANS 
-- relation between Taylor coefficients current and next step 
a0’? = aO + al + a2; 
al’ 
a2’ 


ou 

(| | 

on 
oO 

os. 

* 

A) 

Hare 

I 

tal 

on 

* 

A) 

NO 


-- x(t) is given by a Taylor form of degree 2 
x’? = a0’? + al’*t + a2’*(t°2); 
-- note the implicit existential quantification of t 


-- increment the counter by 1 after each time frame 
counter’ = counter + 1; 


TARGET 

-- state to be reached in bounded time 
x > 1.2 and 
counter <= 9; 


Listing 2: The encoding of Example 1 in iSAT3. 


Example 2 Consider the same DDE equation as Example (1) with the same 
initial condition, but for solving the conjectured safety property of (bounded) 
until operator (x < 1.2) Uoig) (x < 1.0). 
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This time, we employ four Boolean helper variables, of which iSAT’s 
BMC mode will instantiate a fresh copy in each step: 


1. Boolean state variable b records a sufficient condition for « < 1.2 being 
true throughout the current time frame in the sense that b is true only 
if x(t) < 1.2 holds for each time instant t in the current time frame 
(cf. lines 28 and 50 in Listing 3); 


2. Boolean state variable c records a sufficient condition for x < 1.0 being 
true throughout the current time frame (cf. lines 31 and 52); 


3. the Boolean state variable u records a sufficient condition for the 
temporal property (2 < 1.2) Uojio—n) («© < 1.0) being true in the 
current step, with n being the number of the current step (cf. line 54); 


4. Boolean state variable done is a helper variable necessitated by the 
confined expressiveness of the BMC mode, which permits reference 
to current and next states only. It records whether the termination 
condition x < 1.0 has already been true in the past (lines 34 and 55). 


DECL 
-- the range of each variable has to be bounded 
float [-1000, 1000] a0, al, a2, x1, x2; 
float [0,1] ti, t2, xi; 
-- each of the atomic subformulae needs its own 
-- fresh copy of the state variable x and the time 
-- instant t due to the quantifier elimination 


-- define counter for bounded verification problem 
int [0,9] counter; 


-- define Boolean helper variables 

boole b, c, u, done; 
-- b records sufficient condition for x <= 1.2 
-- c records sufficient condition for x <= 1.0 
-- u records sufficient condition for until 


-- done records whether c has been true in the past 


INIT 
xl = 1; 
x2 -= 4: 
-- initialize Taylor coefficients 
aO = 1; 
ail = 0; 


25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
Al 
42 
43 
44 
45 
46 
AT 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
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--initialize b, the sufficient condition of everywhere x <= 1.2 
(not b) -> (x1 > 1.2); 


--initialize c, the sufficient condition of everywhere x <= 1.0 
(not c) -> (x2 > 1.0); 


-- initialize done, the variable memoizing x <= 1.0 
done <-> cs 


-- counter observes the time interval 
counter = 9; 


TRANS 
-- description of the transition system of DDE model 


a0’ = aO + al + a2; 
ai’ = -a0; 
a2’? = -0.5*al - xi *a2; 


-- tracing the bounded until 

-- find witness points 
x1’? = a0’ + al’*ti + a2’*(t1i72); 
x2’? = a0’? + al’*t2 + a2Q’*(t272); 

-- b is sufficient condition for x <= 1.2 throughout time frame 
(not b’?) -> (xi’ > 1.2); 

-- c is sufficient condition for x <= 1.0 throughout time frame 
(not c’) -> (x2’ > 1.0); 

-- recurrence rules for until 
u <-> done or (b and wu’); 


done’ <-> (c’ and counter > 0) or done; -- remembers c 
(counter > 0) -> (counter’ = counter-1); 
(counter = 0) -> (counter’ = 0); 

TARGET 


-- for constructing a counterexample, the until formula ought 
-- to be violated in the initial time instant 
(not u) and (counter = 9); 


Listing 3: The encoding of Example 2 in iSAT3. 


Checking above example for an appropriate unwinding depth of at least 
9, iSAT will report unsatisfiable, which approves absence of a counterexample 
and thus proves the property to be satisfied. 


Example 3 This example (taken from [44]) is an adaptation of Gustafson’s 
model of nutrient flow in an aquarium [17, p. 589f]. It deals with using 


ao 


FOOMOAN DOB WNHH 


100 P. Nazier Mosaad, M. Franzle, B. Xue 


a radioactive tracer for the food chain consisting of two aquatic plankton 
varieties drifting with the currents. The variables in this three-dimensional 
system reflect the isotope concentrations in the water, a phytoplankton species, 
and a zooplankton species, respectively. The original model was an ODE 
model; a concise model would presumably have to use PDE (partial differential 
equations) to model spacial variations and the necessary drifts of species in 
the predator-prey part of the food chain; our DDE model here is a compromise 
between these two extremes. Therefore consider the three-dimensional linear 
DDE 


“3 6. 3B ; 
1 6 —5 


with initial condition £({0, 1]) = [10,0,0] and a conjectured MITL formula 
specifying the distance between the isotope concentrations of two aquatic 
plankton varieties always stays below 10 in a bounded time [0,50], 2.e., 
[0,50] | Le LS |< 10. 


Using Taylor models of degree 1, we calculate the operator relating 
successive parameter vectors to be 


i oe 0 oOo oO 
He —34 6 BE Be BR 
_{o0 0 [gee 0 0 
PII Dee GH ey, 18) . nO He 
0 oO oO oe a oe 
i 3 6 63 —5 —5€3 


In this example, the solver outputs that the system is safe, which means 
that any state satisfying the target property is unreachable within depth 50 
w.r.t. the over-approximation model of the DDE. This constitutes a rigorous 
proof that the system actually satisfies the given property. 


DECL 

-- the range of each variable has to be bounded 
float [-1000, 1000] a01, aii, a02, a12, a03, a1i3, x1, x2, x3; 
float [0,1/100] t, xi1t, xi2, xi3; 


-- define counter for the bounded verification problem 
int [0,49] counter; 


INIT 
-- initial values for the three components of the state 
x1 = 10; 
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x2 = 0; 
x3 = 0; 
-- initialize Taylor coefficients 
a01 = 10; 
alt = 0; 
a02 = 0; 
al2 = 0; 
a03 = 0; 
a1l3 = 0; 
-- initialize the counter 
counter = 0; 
TRANS 
--description of the transition system of DDE model 
x1’ = aOi’ + ali’*t; 
x2’ = a02’ + al2’*t; 
x3’ = a03’ + al3’*t; 
aQ1’? = aO1 + (1/100)*a11; 
aii’? = ((-3)*a01) - ((3*xi1l)*a11) + (6*a02) 
+ ((6*xil)*a12) + (5*a03) + ((5*xil)*a13); 
a02’ = a02 + (1/100)*a12; 
al2’? = (2*a01) + ((2*xi2)*a11) - (12*a02) - ((12*xi2)¥*a12); 
a03’ = a03 + (1/100)*a13; 
a13’ = aO1 + (xi3*aii) + (6*a02) + ((6*xi3)*a1i2) 
- (5*a03) - ((5*xi3)*a13); 
-- increment the counter by 1 for each time frame 
counter’ = counter + 1; 
TARGET 
-- state to be reached in bounded time 
abs(x2-x3) > 10 and counter <= 49; 


Listing 4: The encoding of Example 3 in iSAT3. 


Along these lines, we should point out that the above verification proce- 
dure for temporal specifications could obviously fail, in the sense of providing 
false negatives and corresponding counter-examples, due to excessive over- 
approximation of the DDE’s solution or due to the conditions used for MITL 
satisfaction not being necessary, but just sufficient. As the former problem 
would be induced by selecting an insufficient bound on the degree of the 
Taylor forms, one could simply select a higher degree. Therefore it should, 
however, be clear that the negative verdict actually is spurious due excessive 
over-approximation. Automatic methods to check whether the reported 
counterexample is spurious or not remain to be developed. One solution to 
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this problem could be by using counter-example guided abstraction refine- 
ment (CEGAR) [8] for enhancing the over-approximation model. Another 
solution could be by refining the over-approximation Taylor model using 
sensitivity analysis, and hence enhancing the constructed conditions on the 
model. That latter solution aims at eliminating the wrapping effect due to 
the dependency issue in interval arithmetic. A pertinent algorithm based on 
this idea is currently under development and will be exposed in future work. 


6 Conclusion and Future Work 


In this paper, we have elaborated a method to verify/falsify temporal specifi- 
cations of time-delay systems modeled by a simple class of delay differential 
equations (DDEs) with a single constant delay. Several dynamical systems 
can be modeled by DDEs with a single constant delay as in biology [16, 26], 
optics [19], economics [41, 42], ecology [12], to name just a few. As requi- 
rements specification language, we have exploited metric interval temporal 
logic (MITL) [1] with continuous-time semantics on the solutions of the 
DDEs. We have built our method around a fixed degree interval-based 
Taylor over-approximation technique [44] in order to provide a safe enclosure 
method for DDEs, thereby obtaining timed state sequences spanned by the 
piecewise valid Taylor coefficients. In this way, the continuous semantics of 
the MITL formulae is reduced to a time-discrete problem on timed state 
sequences in terms of Taylor coefficients. Then, we have devised sufficient 
conditions on these timed state sequences recovering the continuous-time 
interpretation of MITL on the actual solutions of the DDEs. To achieve 
this, we have first built sufficient conditions for validation of the atomic 
predicates over time frames of the Taylor over-approximation model of DDE. 
We have then extended this approach to arbitrary bounded MITL formulae 
in negation normal form. Exploiting this as a tableaux or using a related 
encoding as a bounded model checking (BMC) problem, we could employ 
arithmetic SMT solver addressing (a.o.) polynomial arithmetic as a tool 
able to automatically provide certificates of temporal properties for DDEs. 
In our case, we have used the iSAT3 solver, which is the third implementa- 
tion of the iSAT algorithm [13]. In very first experiments on simple DDEs, 
the iSAT3 solver proved able to solve the temporal properties expressed in 
MITL formulae, thereby safely determining satisfaction of the formulae in 
an over-approximation setting. We were able to verify formulae of temporal 
logic also involving Boolean connectives and temporal modalities, like the 
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(bounded) until operator. 

We have presented some examples to demonstrate our method. The 
soundness of the method is guaranteed due to the over-approximation em- 
ployed in DDE enclosure by Taylor forms and the sufficient conditions of 
determining the truth values of the atomic propositions over the time frames. 
Such over-approximation may, however, provide spurious counterexamples 
in case of a failing verification attempt, which ought to be disambiguated 
from true counterexamples. To resolve that ambiguity in case of a negative 
verdict, as a future work, further techniques remain to be developed. We may 
build our idea on the general counter-example guided abstraction refinement 
(CEGAR) technique [8]. Another solution could be by refinement based on 
sensitivity analysis. A pertinent algorithm based on sensitivity analysis is 
currently under development and will be exposed in future work. 


In control applications, one may furthermore want to combine delayed 
feedback, as imposed a.o. by networked control, with immediate state feed- 
back modeled by ordinary differential equations (ODEs). We may investigate 
in more detail some algorithms to handle such cases in the near future. The 
main idea is based on a layered combination of Taylor-model computation 
for ODE, e.g., [31], with the ideas imposed in [44] for DDE. In this way, we 
may extend our method exposed herein to verify the temporal properties 
of dynamical systems modeled by the combination of ODE and DDE. In 
subsequent steps, we plan to extend the method even further to more general 
kinds of DDE, like DDE with multiple different discrete delays, DDE with 
randomly distributed delay, or DDE with time-dependent or more generally 
state-dependent delay [23]. Finally, we would like to point out that in this 
paper, we essentially have presented a verification method based on model 
checking to design time-delay continuous systems modeled by a simple class 
of DDEs. This method may also be used in interactive proofs and stepwise 
refinement of hybrid systems featuring delayed feedback, akin to the methods 
developed for traditional hybrid systems [2, 5]. 
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